Comparison of the toy_06 model

In [1]:
import os, itertools
import pandas as pd
import numpy as np

from bokeh.plotting import figure, output_file, show
from bokeh.io import output_notebook
from bokeh.palettes import Dark2_5, colorblind
from glob import glob

from bokeh.layouts import layout, column, row
from bokeh.models.widgets import Tabs, Panel
from bokeh.io import curdoc
from bokeh.plotting import figure
from scipy import interpolate, integrate

from IPython.core.display import display, HTML
display(HTML("<style>.container { width:100% !important; }</style>"))
In [2]:
def read_magnitudes(fname):
    COL_NAMES = ['time', 'U', 'B', 'V', 'R', 'I', 'J', 'H', 'K']
    mag_table = pd.read_csv(fname, delim_whitespace=True, comment='#', names=COL_NAMES)
    mag_table.iloc[:, 1:] = mag_table.iloc[:,1:].replace(99.999, np.nan)
    return mag_table
    
def read_spectrum(fname):
    col_names = ['wavelength']
    with open(fname, encoding='utf8') as fh:
        for i, line in enumerate(fh):
            if i==2:
                times = list(map(np.float64, line.strip().split(':')[1].split()))
    
    spectrum = pd.read_csv(fname, delim_whitespace=True, comment='#', header=None, names=col_names + times)
    return spectrum

def read_tgas(fname):
    col_names = ['velocity']
    with open(fname, encoding='utf8') as fh:
        for i, line in enumerate(fh):
            if i==2:
                times = list(map(np.float64, line.strip().split(':')[1].split()))
    
    spectrum = pd.read_csv(fname, delim_whitespace=True, comment='#', header=None, names=col_names + times)
    return spectrum
In [3]:
# uncomment the following line to download the code comparison stuff
#!wget http://opensupernova.org/~wkerzend/files/ccompare1_models.tgz
#!tar zxf ccompare1_models.tgz

Run interactively on Google Colab on the following link

https://colab.research.google.com/github/sn-rad-trans/code-comparison1/blob/master/plotting_toy01.ipynb

In [4]:
CODE_NAMES = ['artis', 'cmfgen', 'sedona', 'stella', 'sumo', 'supernu', 'urilight']
CODE_COLORS = {name:color for name, color in zip(CODE_NAMES, colorblind['Colorblind'][8])}

Comparison of pseudo bolometric luminosity

Integration of the synthetic spectra excluding gamma rays ($>50 Angstrom$)

In [5]:
output_notebook()


# create a new plot with a title and axis labels
fig_lbol  = figure(title="Comparison of Lbol", x_axis_label='time [days]', y_axis_label='Log Lbol')
fig_edep  = figure(title="Comparison of energy deposition", x_axis_label='time [days]', y_axis_label='Log Edeptot')
fig_econ  = figure(title="Energy conservation", x_axis_label='time [days]', y_axis_label='(int t * lbol * dt)/(int t * edep * dt)')

lbol_data = {}
for code_name in CODE_NAMES:
    if code_name == 'sedona':
        fname = 'toy01/lbol_edep_toy01_{0}.txt'.format('sedonaexpansion')
    else:
        fname = 'toy01/lbol_edep_toy01_{0}.txt'.format(code_name)
    if not os.path.exists(fname):
        continue
    lbol = pd.read_csv(fname, delim_whitespace=True, comment='#')
    lbol.columns = ['time', 'lbol', 'edep']
    lbol['log_lbol'] = np.log10(lbol.lbol)
    lbol['log_edep'] = np.log10(lbol.edep)
    lbol['econ']= np.hstack(([np.nan], integrate.cumtrapz(lbol.lbol*lbol.time, lbol.time) / integrate.cumtrapz(lbol.edep*lbol.time, lbol.time)))

    fig_lbol.line(lbol.time, lbol.log_lbol, legend=code_name, line_width=2, color=CODE_COLORS[code_name])
    fig_edep.line(lbol.time, lbol.log_edep, legend=code_name, line_width=2, color=CODE_COLORS[code_name])
    fig_econ.line(lbol.time, lbol.econ, legend=code_name, line_width=2, color=CODE_COLORS[code_name])
    lbol_data[code_name] = lbol

fig_lbol.legend.location = "top_right"
fig_lbol.legend.click_policy="hide"
fig_edep.legend.location = "top_right"
fig_edep.legend.click_policy="hide"
fig_econ.legend.location = "top_right"
fig_econ.legend.click_policy="hide"


# show the results
show(column([fig_lbol, fig_edep, fig_econ]))#
Loading BokehJS ... 
/Users/wkerzend/miniconda/envs/sn-rad-trans/lib/python3.7/site-packages/pandas/core/series.py:853: RuntimeWarning: divide by zero encountered in log10
  result = getattr(ufunc, method)(*inputs, **kwargs)

Comparing gas temperatures

In [6]:
## Reading the data
tgas_data = {}
for code_name in CODE_NAMES:
    if code_name == 'sedona':
        fname = 'toy01/tgas_toy01_{0}.txt'.format('sedonaexpansion')
    else:
        fname = 'toy01/tgas_toy01_{0}.txt'.format(code_name)
    if not os.path.exists(fname):
        continue
    try:
        tgas_data[code_name] = read_tgas(fname)
    except UnicodeDecodeError:
        pass
    except IndexError:
        pass
In [7]:
TGAS_TIMES = [5, 10, 15, 20, 30, 50, 100, 200, 300, 400]
output_notebook()
interpolators = {}
for name in CODE_NAMES:
    if name not in tgas_data:
        continue
    interpolators[name] = interpolate.interp1d(tgas_data[name].columns[1:], 
                                               tgas_data[name].iloc[:, 1:].values, 
                                               fill_value=np.nan, bounds_error=False)
all_tabs = []
for time in TGAS_TIMES:
    fig = figure(x_axis_label='velocity [km/s]', y_axis_label='Tgas [K]')
    for name, interpolator in interpolators.items():
        velocity = tgas_data[name].velocity
        fig.line(velocity, interpolator(time), legend=name, line_width=2, color=CODE_COLORS[name])
    fig.legend.location = "top_right"
    fig.legend.click_policy="hide"
    lout = layout([[fig]])
    tab = Panel(child=lout,title='t={0:.0f}'.format(time))
    all_tabs.append(tab)
        
show(Tabs(tabs=all_tabs))
Loading BokehJS ...

Comparing the magnitudes

In [8]:
## Reading the data
mag_data = {}
for code_name in CODE_NAMES:
    if code_name == 'sedona':
        fname = 'toy01/mags_toy01_{0}.txt'.format('sedonaexpansion')
    else:
        fname = 'toy01/mags_toy01_{0}.txt'.format(code_name)
    if not os.path.exists(fname):
        continue
    try:
        mag_data[code_name] = read_magnitudes(fname)
    except UnicodeDecodeError:
        pass
    except IndexError:
        pass
In [9]:
output_notebook()

MAG_NAMES = list(mag_data.values())[0].columns[1:]
all_tabs = []
for mag in MAG_NAMES:
    fig = figure()
    for name, mag_table in mag_data.items():
        fig.line(mag_table.time, mag_table[mag], legend=name, line_width=2, color=CODE_COLORS[name])
        fig.y_range.flipped = True
    fig.legend.location = "top_right"
    fig.legend.click_policy="hide"
    lout = layout([[fig]])
    tab = Panel(child=lout,title=mag)
    all_tabs.append(tab)
        
show(Tabs(tabs=all_tabs))
Loading BokehJS ...

Comparing the colors

In [10]:
output_notebook()

MAG_NAMES = list(mag_data.values())[0].columns[1:]
all_tabs = []
for mag1, mag2 in zip(MAG_NAMES[:-1], MAG_NAMES[1:]):
    fig = figure()
    for name, mag_table in mag_data.items():
        fig.line(mag_table.time, mag_table[mag1] - mag_table[mag2], legend=name, line_width=2, color=CODE_COLORS[name])
        fig.y_range.flipped = True
    fig.legend.location = "top_right"
    fig.legend.click_policy="hide"
    lout = layout([[fig]])
    tab = Panel(child=lout,title='{0}-{1}'.format(mag1, mag2))
    all_tabs.append(tab)
        
show(Tabs(tabs=all_tabs))
Loading BokehJS ...

Comparing the spectra

In [11]:
## Reading the data
spectral_data = {}
for code_name in CODE_NAMES:
    if code_name == 'sedona':
        fname = 'toy01/spectra_toy01_{0}.txt'.format('sedonaexpansion')
    else:
        fname = 'toy01/spectra_toy01_{0}.txt'.format(code_name)
    if not os.path.exists(fname):
        continue
    try:
        spectral_data[code_name] = read_spectrum(fname)
    except UnicodeDecodeError:
        pass
    except IndexError:
        pass
In [12]:
SPECTRAL_TIMES = [5, 10, 15, 20, 30, 50, 100, 200, 300, 400]
output_notebook()
interpolators = {}
for name in CODE_NAMES:
    if name not in spectral_data:
        continue
    interpolators[name] = interpolate.interp1d(spectral_data[name].columns[1:], 
                                               spectral_data[name].iloc[:, 1:].values, 
                                               fill_value=np.nan, bounds_error=False)

all_tabs = []
for time in SPECTRAL_TIMES:
    fig = figure()
    for name, interpolator in interpolators.items():
        wavelength = spectral_data[name].wavelength
        fig.line(wavelength, interpolator(time), legend=name, line_width=2, color=CODE_COLORS[name])
    fig.legend.location = "top_right"
    fig.legend.click_policy="hide"
    lout = layout([[fig]])
    tab = Panel(child=lout,title='t={0:.0f}'.format(time))
    all_tabs.append(tab)
        
show(Tabs(tabs=all_tabs))
Loading BokehJS ...
In [ ]: